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The history force is one of the hydrodynamic forces which act on a particle moving 
through a fluid. It is an integral over the full time history of the particle’s motion and 
significantly complicates the equations of motion (accordingly it is often neglected). We 
present here a study of the influence of this force on particles moving in a turbulent 
flow, for a wide range of particle parameters. It is shown that the magnitude of history 
force can be significant and that it can have a considerable effect on the particles’ slip 
velocity, acceleration, preferential concentration and collision rate. We also investigate 
the parameter dependence of the strength of these effects. 


1. Introduction 


The advection of finite-size particles by a fluid flow plays an important role in a large 
number of natural and industrial situations, a few examples are: cloud microphysics, 
combustion processes and formation of marine aggregates or planets. The underlying 
flow is turbulent in many cases. When the inertia of the particles is small enough they 
can be accurately modelled as tracers, i.e. particles which follow the fluid flow exactly. 
When this is not the case the particles can deviate from the flow due to their inertia and 
one then speaks of inertial particles. A fundamental difference between these two types 
of particles is that the dynamics of tracers preserves phase space volume (when the flow 
is incompressible) whereas the dynamics of inertial particles is dissipative, thus allowing 
for the existence of attractors. As a direct consequence, inertial particles can accumulate 
in certain regions of the flow - an effect termed preferential concentration - while tracers 
stay homogeneously distributed for all times. 

The equations of motion for a particle in a fluid flow can (under certain assumptions) 
be systeniaticalIv^erivedfron^the_^Iavier^Stokes,^yuatioiij,^iSjias_^eei^^niie_^^Maxey 


& Riley HQS,*!! ): Clatismoll " IQS.jl l (see the references therein for a historical overview of 
the literature before 1983). These equations of motion are comprised of different terms 
which makes it possible to decompose the total hydrodynamical force on the particle 
into different contributions like the Stokes drag, the added mass effect and the history 
force. The latter is an int egral over th e whol e history of the particle’s motion - a memory 
effect. Boussinesa ( 18851 ) and [Basset ( 18881) were the first to point out the existence of 
the history force, which is therefore often called the Boussinesq-Basset force. 

Before we proceed to the full equation, let us present an intuitive argument for the 
existence of the history force and try to guess its form. Consider a spherical particle 
moving through a still fluid with a constant velocity. If the Reynolds number is small, we 
can use the linearised Navier-Stokes equation and the flow is then given by the Stokes 
flow around a sphere. Now imagine a sudden change dup in the velocity of the particle 
at time t = 0. This disturbance will lead to the generation of vorticity at the particle’s 
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surface, which will decay diffusively and thus as l/v^ in the vicinity of the particle. 
We might expect that the force on the particle due to this decaying part of the flow 
will be proportional to dvo/v^- This is already a memory effect because the particle 
“remembers” the disturbance at t = 0 in all following instants. The particle obviously 
does not have a memory itself but rather the information is stored in the flow around it. 
Now imagine a general motion of the particle as a series of velocity jumps dt;(r). Then 
the superposition of the flows due to these disturbances will produce a force proportional 
to 




1 dv 

y/t-T dr 


dr 


\/t — T '' ’ Jq \/t — t dr ^ ^ 

(a superposition is admissible because of the small Reynolds number and thus the ap¬ 
proximate linearity of the Navier-Stokes equation). Up to a constant prefactor (If.Ill is the 
history force on a particle moving in a still fluid. It describes the effect of the decaying 
parts of the flow which were generated by the acceleration of the particle at previous 
instants (on the particle at the present instant). 

In the case of a general fluid flow the motion of a spherical particle is described by the 
following equation 


dv 

dt 


/oDm 1 / ^ 



(1.2) 


This is the version derived bv iMaxev fc Ril^ ( 1983 1: Gatignoj (Il983h wit h the slightly 
different (and widely used) form of the added mass term by Auton et all (1981). Here 
V = dx/dt, u = u(x,t) and ^ = dtu + u ■ Vm denote the particle velocity, the fluid 
velocity and the fluid acceleration at the position of the particle; u is the undisturbed 
flow, without the particle’s presence. There are two parameters in (11.211 . the first is the 
density parameter 


where g = Qp/gf is the particle’s density gp normalized by that of the fluid g /. When 
g > 1 the particle is heavier than the fluid and 0 < /3 < 1; when < 1 it is lighter than 
the fluid and 1 < /3 < 3. The second parameter is the particle response time 


1 

3/3 V 


(1.4) 


where r is the radius of the particle and v the kinematic viscosity of the fluid. We have 
omitted here the so-called Faxen corrections and the influence of gravity. The terms 
appearing on the right-hand side of (HU are the pressure gradient (which contains a 
contribution from the added mass term), the Stokes drag and the history force. The 
latter is an integral over the whole history of the particle and describes the effect of the 
decaying disturbance flow generated by the particle at earlier times on the particle at 
the present time, in harmony with the intuitive argument described above. The history 
force is a viscous effect and is sometimes referred to as the unsteady drag. 

In turbulence the relevant time scale to which the particle response time Tp should be 
compared is the Kolmogorov-time = \Jvje, where e is the mean energy dissipation. 


f There is actually another effect of the disturbance duo: the added mass effect, which is 
proportional to duo/dt and describes the immediate pressure response of the fluid. We neglect 
it here for the sake of a concise illustration of the history force. 
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The Stokes number 

St='^ (1.5) 

Trj 

is the ratio of these two characteristic time scales and is a measure for the importance 
of the particle’s inertia. Throughout the paper we will use the Kolmogorov scales of 
turbulence rj = (length), (time), = r]/Tr^ (velocity) and a,, = Ur^jTr^ 

(acceleration) to normalize dimensional quantities. These scales are the smallest scales 
of the turbulent motion and thus relevant for small particles, which “live” on the smallest 
scales, e.g. for the particles considered here Tp/r^ and r/rj are at most on the order of a 
few times unity. 

Note that in (II.2|) the form of the history force is non-standard: the derivative is outside 
the integral. This form is actually more general as it is valid for initial conditions with 
v{0) ^ tt(0). The following relation (derived via partial integration) shows the connection 
between the generalized form and the standard form 


f f ^ dr + 

dt Jo y/t-T Jo d/t-T d/i 



( 1 . 6 ) 


Th e term (t>(0) — u(0 ) ) / y/i should b e added in m when i;(0) ^ 'ir(O) as pointed out 
by iMichaelides ( 1992ll : iMaxevI ( 1993l l. Thus the left-hand side in (II.6p is a form of the 
history force which is valid for any type of initial condition. This generalized form is also 
exactly t he definition of th e fractional derivative of Riemann-Liouville type of order 1/2 
(see, e.g. IPodlubnvi (jl998)). Th i s conn ection to fractional derivatives and relation (jl.6l) 
was first pointed out by Tatom ( 1988l l. We will use the fractional derivative notation to 
abbreviate the history force in the following. 

The history force has been frequently neglected in applications of the equation of 
motion 1121). This is often done to simplify the problem: the history force turns the 
equation of motion into an integro-differential equation and thus makes it much more 
difficult to reason about the solutions. For example, the existence, uniqueness, regularity 
and asymnotiWnsofthe_sohitions£yT2JlJiave^eei^rail^^ecenK^,studied^^^Farazmand 
& Haller ( 2015l l: lLanglois et all (120141 ). Another difficulty is the computation of numerical 
solutions with the history force: on the one hand there is the singularity of the integrand of 
the history force which impedes an accurate numerical approximation and on the other 
hand there is the necessity to recompute the history integral for every new time step 
which causes high numerical costs. The first problem can be resolved by an appropriate 
treatmen t of th e history force and, by now, higher-order integration schemes are available 
( Daitchd ( 20131 ) ). The second problem is inherent to the dynamics wi th memory and can 


be att enuated only by an approximation of the history kernel (see e.g. Ivan Hinsberg et al. 

(mil)). 


Druzhinin & Ostrovskv 

(I 1994 I); Coimbra & Rangel (' 

1998) 

Hill (2005 

1: Candelier et al. 

(20041. 20141): Lim et al. 

(I 2 OI 4 I) and experimentallv bv 

Mordant & Pinton] (2000|)^ Abbad 

& Souhar ( 20041) ;ICoimbra et al\ ( 20041) ; iToegel et al\ 

200f^):ICarbin et a,l\ (2009^. It was 


shown that memory can be quite important, for example the experimental studies showed 
that the history force is necessary (in these cases) for a match bet ween experiment and 
theory . Chaotic dynam i cs of i nertial particle s have b een studied by Yannacopoulos et al 


tneory . t^naotic dynam i cs ot i nertial particle s nave b een studied by Yannacopoulos et al. 
( 19971 ): Daitche fc Tel ( 201 il l: Guseva et al\ ( 2013h : IPaitche k. Tel ( 20l4l showing that 


the history force can qualitatively change the dynamics, for example, the history force 
reduces the tendency for accumulation and can change the n a ture of attractors from 
non-chaotic to chaotic and vice versa. iReeks fc McKeel ( 1984 ): Mei et al. 1 199lh : van 
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Aartrijk & Cl ercx (1201(111 have shown t hat memory can affect the dispersion of particles 
in turbulence. ICalzavarini et al. ( 2012 1 studied the influence of Faxen corrections, non¬ 
linear drag and history force on the statistical properties of large neutrally buoyant 
particles {r/rj G [3, 30]) moving in a turbulent flow and found a significant influence 
of all three terms (especially for particles considerably larger than t]). The importance 
of differen^_forces_^cting_^^aDarticle_ia^turbulence_ias_^eea^studie^^^Armenio & 
Fiorotto ( 200lh : lOlivieri et all (l2014h showing that the magnitude of the history force 
can be significant and that it can also alter the contrib ution o f the other forces. Olivieri 


et al. ( 2014l l also found that memory can reduce preferential concentra tion. We will 


present a detailed comparisons with the results by lOlivieri et all (1201J) in section [51 


where we find a considerable difference to our results. 

One of the conditions for the derivation of (11.211 is that the 


Rep = r\v — u\ /v is small. Lovalenti fc Bradvl (|l993 1 and iMeil ( 199-dl have de veloped 


ar ticle R eynolds number 


extended versions of (II.2|) for finite Rcp (up to a few hundreds bv iMeil (|l994ll l. Both 
versions contain significantly more compl icated forms of the history force; Mei’s version 
also includes a nonlinear form of the drag. lMaxev et aJ. 1 1996h found that up to Rcp Ri 17 
the standard form (11.21) “may be quite adequate in practic e even though not justified by 
theory”. In an experimental studv iAbbad fc Souhan ( 2004 1 considered particle Reynolds 
numbers of up to 0.5 and found the standard form to remain valid. In our simulations 
the typical Rep is smaller or on the order of unity (except for p = 0 and St > 2, when it 
becomes larger than 3). We thus assume that the standard version (11.21) is applicable in 
our case. As we want to concentrate on the effects of the history force we do not consider 
here any further corrections, like e.g. the lift force. In a comparison of the magnitudes of 
the Faxen corrections and the history force we found that the former are much smaller 
in most cases considered here, see appendix [^ Thus we neglect the Faxen corrections in 
the following. 

The aim of the present paper is to provide a comprehensive study of the role of the 
history force for the motion of particles in a turbulent flow. To this end we prepared 
numerical simulations for a wide range of particle parameters, which will be investigated 
for effects of memory on different statistical particle properties. An important question we 
want to answer here is for which particle parameters the effects of memory are important. 

The paper is structured as follows: The next section will present details on the numer¬ 
ical simulations. The magnitude of the history force will be studied in section El followed 
by investigations of its effect on the slip velocity v — u (section [4]), the acceleration (sec- 
tionE]), preferential concentration (section al) and collision rates (sectionO. Section|5]will 
detail a comparison with the recent work bv lOlivieri et al\ ( 2014h . followed by a summary 
and discussion in section [9] Appendix [^ presents a comparison of the history force and 
the Faxen corrections. A note on the case of neutrally buoyant particles, which is not 
considered in the main part of text, is given in appendix |B] The numerical scheme for 
the integration of particle trajectories with the history force and the method of forcing 
the turbulence are described in the appendices |C] and |D1 respectively. 


2. Numerical simulations 

The turbulent flow is generated in a triply periodic box by a large scale forcing (as 
described in appendix IQ. We solve the vorticity equation, which is equivalent to the 
incompr essible Naiver- S tokes e quation, b y a sta ndard dealiased Fourier-pseudo-spectral 
method ( Canuto et al. (Il987h: IHou fc Lil (|2007l ll with a third-order Runge-Kutta time¬ 
stepping scheme ( Shu fc Osheil ( 1988h l. The particle trajectories are integrated with a 
specialized scheme which treats the history force appropriately; it is a modified version of 
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Rcx L / Tj ^/rj j Tj Tsijrj^^Trj TjrT^ ^ijr-q UYmsIUjj hi 

113 632 157 20.9 1.23 997 29.1 0.0199 5.40 512® 

Table 1. Parameters of the simulated turbulent flow: Taylor Reynolds number Re\ = Xu^ms/n, 
size of the periodic box Lbox, integral scale L = u®nia/e, Taylor microscale A = Wrms-\/l5z.^/e, size 
of a grid cell A®, length of the whole simulation Tsim, large-eddy turnover time T = L/uims, 
time step At, root-mean-square of the velocity Urms = x/{iiP)/3, number of grid points A®. All 
dimensional quantities are given in multiples of the corresponding Kolmogorov scales. 


the third-order scheme developed bv IPaitche ( 2013 1 and is described in appendix [Cl The 
values of Euleriaii quantities, which are present on a grid, are obtained at the particle 
positions through tricubic interpolation. The turbulent flow is statistically homogeneous, 
isotropic and stationary, with the flow characteristics depicted in table HJ 

For this study a number of simulations with different particle parameters have been 
prepared. For the density g the values 1000, 100, 10, 2, 0.5, 0 were chosen, g = 1000 
corresponds to water droplets moving in air, relevant for, e.g., cloud microphysics. Density 
ratios 10 and 2 are close to those of metals and sand in water, p = 1 is the case of 
neutrally buoyant particles, which is not discussed in the main text but in appendix |BJ 
g = 0 corresponds to air bubbles in water. The Stokes number has been varied in the 
range [0.1, 3.0]. It is typically in this range where interesting properties of inertial particles 
appear. Also, this range is constrained from above by the limitations of (11.21) (the particle 
Reynolds number grows with St) and from below by the time step of the simulation (it 
should be significantly smaller than Tp). For every combination of the values of g and 
St a simulation of particles with and without the history force has been prepared. The 
number of simulated particles for every parameter combination is Np = 10®, except 
for the investigations in sections |6| and |7| where Np = 5 • 10® (in this case the length 
of the simulations is Tsim = TOOr,, = 24T). The initial particle positions have been 
chosen randomly and homogeneously distributed in space; the initial particle velocity 
is that of the fluid at the particle’s position. Statistical quantities (e.g. averages) have 
been obtained by sampling over the particle ensemble and over time. To assure that the 
particles have equilibrated with the flow, this sampling is started after an initial period 
of five large-eddy turnover times. 


3. Forces 

We start by comparing the magnitudes of the different forces acting on a particle. 
In the following, the moduli of the particle acceleration and of the three terms on the 
right-hand side of (|1.2I1 will be denoted by a, ap (pressure gradient), as (Stokes drag) 
and an (history force), respectively. 

Figures [1^ andjUD show the probability density functions (PDFs) of ap, as and ap 
for particles heavier (g = 2) and lighter (g = 0) than the fluid. In both cases ap is the 
dominant contribution, followed by as and ap. The Stokes drag and the history force 
are of similar magnitude showing that the history force is a viscous effect which can be 
as important as the drag. In these particular cases we find (ap) / (as) = 66% in figurejT^ 
and (ap) / (as) = 150% in figure (Ud. Additionally the PDF of the history force has 
longer tails than that of the Stokes drag (see hgure |T^ and also note that the x-axis is 
logarithmic), i.e. the history force has more extreme events. 

Let us now consider the relative importance of the history force and compare it to the 
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ai/{a) 


ai/{a) 


Figure 1. PDFs of the different accelerations at due to the pressure gradient [i = P), the Stokes 
drag (i = S) and the history force [i = H). The particle parameters are (a) q = 2, St = 0.5 
{r/rj = 0.95) and (b) g = 0, St = 0.5 (r/ri — 2.1) . The vertical lines show the positions of 
the averages: (ap) / (a) — 0.62, {as) / (a) = 0.25, (ur) / (a) = 0.17 in (a) and (ap) / (a) — 2.4, 
(as) / (a) = 0.62, (oh) / (a) = 0.93 in (b). 


Stokes drag by means of the ratio oh/ as- Before proceeding to the numerical simulations, 
let us try to obtain an estimate of this ratio. To this end we apply the approximation 


dt J 


1/2 


{v — u) 


a 


(v-u), 


(3.1) 


which is motivated by the fact that the history force is a (fractional) time-derivative. The 
choice of as the characteristic time scale oi v — u seems natural as the particle “lives” 
on the small scales and should follow the flow to some degree (except for very large St). 
The constant a is expected to be on the order of unity. With this approximation and the 
definition of a n and as we obtain a n estimate for the relative magnitude of the history 
force (see also Daitche fc Tell ( 20141 )): 


an 

as 



(3.2) 


Thus we expect the importance of the history force to scale solely with the particle size 
and to be independent of its density. 

Figured^ shows (an) / {as), measured in the turbulence simulations, as a function of 
r/r]. We see a collapse of the data for different densities, which confirms our expectation 
that the particle size is the determining parameter for the relative magnitude of the 
history force. (When looking more closely we actually see a very weak dependence on 
the density, the collapse is not perfect.) A fit of our prediction (|3.2I) to the data yields 
a = 0.69 and we see that this simple estimate holds quite well. 

Figurej^b shows (an) / (as) as a function of the Stokes number. In this case there is no 
collapse for different densities and the relative magnitude of the history force decays with 
growing density. This tendency is true for a fixed Stokes number. Because St depends on 
both the density and the particle size, increasing g and holding St fixed means that we 
implicitly decrease the particle size. Thus the dependence on g in figure [2 ]d is actually 
a disguised dependence on r. That there is no “real” dependence on g can be concluded 
from figure [5^. 
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(a) (b) 




Figure 2. The ratio of the average magnitude of the history force and the Stokes drag as a 

function of (a) r/p and (b) St. 


(a) (b) 



Figure 3. (a) PDFs of the slip velocity for particles with and without memory for 
Q = 2 and St = 0.5 {r/p = 0.95). The vertical lines show the position of the averages: 
{\v - u\ /u^}n,emory = ^.24 and (lu - u\ memory = (b) The average slip velocity 

with (filled symbols) and without memory (unfilled symbols) as a function of the particle size 
and for different particle densities. 

We are now able to make an informed categorization of the importance of the history 
force based on its relative magnitude. For r/p around 10“^ or less its contribution is 
very small (less than 1%) and should be negligible in most cases. For r/p around 0.1 the 
contribution is around 7%; while not negligible the history force is expected to play a 
minor role. For r on the order of p the contribution of the history force is on the order 
of the Stokes drag and is thus expected to be important for the particle motion. 

4. Slip velocity 

Now that we can estimate the magnitude of the history force, let us turn to its effects 
on the particle dynamics. A basic property of inertial particles is that they can deviate 
from the fluid flow (in contrast to tracers). It is thus natural to study the effect of the 
history force on this deviation, namely on the slip velocity v — u. Figure[3K shows the PDF 
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r/rj St 


Figure 4. The ratio of the average slip velocity with and without memory as a function of (a) 

r/ri and (b) St. 


of |p — it| with and without the history force (denoted by “memory” and “no memory” 
respectively). Memory leads to a reduction of the slip velocity; in this case the mean slip 
velocity is reduced by 35% in comparison to the case without memory. Also, the tails of 
the PDF become shorter, i.e., strong detachments from the flow occur less frequently with 
memory. Figure[3jD presents an overview of the full parameter space, showing the average 
slip velocity with and without memory. The natural parameter for the slip velocity is 
the Stokes numbeiQ, however to make the role of the size clearer the data are shown as 
a function of r/rj. We see that the effect of the history force (the difference between the 
case with and without memory) increases with r. 

To obtain a compact representation of this effect of memory let us use the ratio 
1'*^ - “ I memory / 1'*^ “ “I no memory > i-®- ^lip velocity with the history force divided by 
the one without. Before proceeding to the data from the simulations let us first try to 
obtain an analytical estimate. Using the equation of motion (jl.2l) and the approximation 
iH) we can estimate the slip velocity as 


Dm dM\ ( 1 / 3/3 

Bt dt ) “V TpTr, 


-1 


H S(-) 


where the case without memory is contained for a = 0. If we assume that the first factor 
on the right-hand side of (|4.1I) is the same for both cases we obtain 


|m- 


- IJ 

I memory 

I 

I no memory 


1 

l + a^' 

V 


(4.2) 


Parameter a was introduced in section |31 where its value was determined from the nu¬ 
merical simulations: a = 0.69. Thus (14.21) provides a prediction for the slip velocity ratio. 
Figure Hi shows the ratio {\v - / (|m - M|)no memory obtained from the simula¬ 

tions along with this prediction. Although not perfect, this estimate - based on rather 
simple assumptions - fits reasonably well to the simulations (note that the black line in 
figure 11^ is not a fit of (j4.2ll to the data; we used the previously obtained value of a). 
An important observation is that the slip velocity ratio collapses for different densities 


f More precisely, the weakly inertial limit (18.611 suggests that it is proportional to S't |/3 — 1|. 
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Figure 5. PDFs of the magnitude of the acceleration for particles with and without memory 
and tracers. The particle parameters are (a) q = 2, St = 0.5 {r/y = 0.95) and (b) ^ = 0, 
St = 0.5 {r/y = 2.1). The vertical lines show the positions of the averages: = 2.1 and 

(“/«’7>memory = 1-9- («/“»7)„o memory = m (a) and (a/a^>^emory = 2-5, (a/o^)^mLory = 3-0 

in (b). 


when considered as a function of r. Although this collapse is not perfect (there seems to 
be a systematic deviation for St > 1), it is still remarkable that the density plays only a 
minor role also for the effect of memory on the slip velocity. When the slip velocity ratio 
is considered as a function of St, there is no collapse at all, see figure HJd. 

This result suggests again that the natural parameter for the effect of the history force 
is the particle size. Due to its simplicity, formula (14.21) . as well as (13.21) . can be very useful 
for estimates of this effect (when a is not known, one might assume a Ri f). 

The findings of this section support the categorization of the memory effect for different 
particle sizes given at the end of the previous section. For r/y < 10“^ the effect is negli¬ 
gible as the relative reduction of the slip velocity ( l-(|u - / (It- - ^Dno memory) 

is below 1%. For r/y Ri 0.1 the reduction is around 5%, i.e. not completely negligible but 
still small. For r /77 Ri 1 it is around 35% and thus considerable. 


5. Acceleration 

Let us now turn to the influence of the history force on the particle acceleration, a much 
studied quantity in turbulence research. Figures[5^ and[5lD show the PDFs of the modulus 
of the acceleration for particles with and without memory and tracers. For heavy particles 
(figure [5^) the history force increases the acceleration: the average increases slightly and 
the tails more notably (note the logarithmic x-axis). For light particles (figure [^r) the 
effect is opposite, the history force decreases the acceleration. The change of the mean 
acceleration is stronger in this case. We can describe these observations in a unified way: 
the acceleration PDF with memory comes closer to that of tracers (compare figures [5^ 
and [SJd) . 

To obtain an overview of the parameter dependence we again turn to ratios. Figure [ 6 ] 
shows the ratio of the (averaged) magnitudes of acceleration with and without memory. 
We see that the history force generally increases the acceleration of heavy particles and 
decreases the acceleration of light particles, as we found in two particular cases above. 
This effect is noticeable but not very large, except for p = 0 and large particle sizes. Note 
that there is no collapse of the ratios shown in figure [ 6 ] in contrast to the slip velocity 
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Figure 6. The ratio of the averaged acceleration magnitudes with and without memory. 



Figure 7. The distribution of light particles (a) with and (b) without memory. Both plots 
show the positions of 5 • particles contained in the slab [0, Lbox] x [0, Lbox] x [0, Lbox/10] at 
t = 291r^ (this corresponds to 10 large-eddy-turnover times). The particles were homogeneously 
distributed at the beginning of the simulation and have the parameters ^ = 0 and St = 0.5. 

ratio shown in figure 0^ and the force ratio shown in figured^. Although the influence of 
memory generally increases with r also for this observable, there is additionally a clear 
dependence on the density. 


6. Preferential concentration 

A much studied phenomenon concerning inertial particles is preferential concentration: 
the tendency of particles to gather in certain regions of the flow. Figure [7] shows the 
distribution of light particles with and without memory. In both cases we can see an 
inhomogeneous distribution of particles, i.e. preferential concentration. We also clearly 
see that the preferential concentration is weaker with the history force (figures [7^ and 
Eb show the same number of particles). 

To obtain an overview of the parameter dependence of this effect we need to quantify 
preferential concentration. There are several possible approaches; we use the correla- 
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Figure 8. The correlation dimension D 2 of the particle distribution in space as a function of 
St for different densities with and without memory. 


tion dimension D 2 as introduced by Grassberger fc Procaccia ( 1983ll . Let P{R) be the 
probability to find a given particle pair with a distance below R, i.e. 

P{R) = {Q{R-\xi-X 2 \)), (6.1) 

where 0 is the Heaviside function and Xi, X 2 the positions of the two particles. Then 
the correlation dimension is defined through 

P{R) (xR°^ R^ 0. (6.2) 

In studies of preferential concentration the radial distribution function g{R) is frequently 
used. It is closely related to P{R) and £>2 {d is the space dimension): 


9iR) « 


(6.3) 


see e.g. iBec et al\ (120051 ). Thus D 2 can also be seen as a characterization of the radial 
distribution function. 

To determine D 2 we measured P{R) in our simulations and fitted a power law in 
the range R/rj G [l0“^,10]. Figure [8] shows the correlation dimension as a function of 























12 


Anton Daitche 



rh 

Figure 9. The change in the correlation dimension Di due to memory as a function of the 

particle size and density. 


the Stokes number for different densities. Witho ut memory we see the ty pical unimodal 


shape of D 2 {St) found in previous studies, e.g. ICalzavarini et al\ (j2008f l. The minimal 
value of 112 , corresponding to maximal preferential concentration, is obtained around 
St K, 1. Concerning the influence of the history force we hnd: for g = 1000 there is hardly 
any effect; for p = 100 there is a small effect; for p = 2 and = 10 the effect is quite 
noticeable; for light particles (gi = 0 and g = 0.5) the effect is strong, e.g. for p = 0 the 
smallest value of D 2 increases by 1.1 when memory is included - a very considerable 
change for a fractal dimension. It is interesting to note that the position of the smallest 
D 2 is shifted by the history force for p = 2; for p = 0.5 and g = 0 there does not seem 
to be a clear minimum with memory but rather an extended plateau. Also note that 
the influence of memory is strongest where D 2 attains its minimum. To sum up, if we 
fix the Stokes number and decrease g, the effect of the history force on D 2 and thus on 
preferential concentration increases and becomes very strong for small densities. 


Note that for a fixed Stokes number we cannot discern the dependencies of this effect 
on r and g. As St depends on both parameters, varying g implicitly varies r; in particular, 
decreasing g increases r, see (IT3])-(fT3D. Thus the dependence on g in hgure [5] might be 
a disguised dependence on r (as is the case in hgure [2]3 for the relative magnitude of the 
history force and in hgure IHd for the slip velocity ratio). Thus we have to consider this 
effect as a function of the particle size and density to discern the dependencies. Figure |9] 
shows the change of the correlation dimension 192, memory ~ 192, no memory as a function 
of these two variables. We see a dependence on both r and g. This is in contrast to the 
behaviour of the relative magnitude of the history force (section [3]) and the slip velocity 
ratio (section |4]), where the determining parameter was solely r and similar to the case 
of the acceleration (previous section). Note also that the dependence on r and g is quite 
non-trivial, e.g. it is non-monotonic in both variables. 

A good way to summarize these findings is as follows: For a fixed Stokes number the 
influence of the history force increases monotonically when the particle density decreases 
or (equivalently) when the particle size increases. However, one has to keep in mind that 
this dependency cannot be attributed to only one of the two parameters. 
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7 . Collision rates 

The collision rate of particles in turbulent flows is important for the understanding of, 
e.g., cloud microphysics or aggregation processes; it is also closely related to preferential 
concentration. One of the main findings in this area is that turbulence can enhance the 
collision rate enormously. We will study here the effect of the history force on the collision 
rate. 

Let Nc{t, At) be the number of particle collisions during a small time interval [t, t + At], 
then the collision rate per unit volume is 


R{t) = 


Ncit, At) 
^LxAr ■ 


(7.1) 


This quantity is independent of t in stationary turbulence, but grows quadratically with 
the particle number density n = Np/L^^^. A quantity where this dependence is factored 
out is the collision kernel 


r = 


2R 


(7.2) 


which we will use in the following. 

To obtain Nc{t, At) we have simulated an ensemble of Np = 5 • 10^ particles for every 
parameter value and saved snapshots of the particle positions every At = 0.2t^ time 
units. The collisions have been determined using the so-called ghost collision approx¬ 
imation, i.e. in our simulations the particles do not interact on contact but simply slip 
through each other. A collision is considered to occur when the distance between two 
particles becomes less then 2r. To determine the occurrence of such events, the particle 
trajectories have been interpolated linearly in between the saved snapshots. 

Figure (Tn] shows the collision kernel as a function of the Stokes number for different 
densities. We see that the inclusion of the history force leads to a decrease in the collision 
rate. This effect is very strong for small densities and decreases for denser particles (when 
St is fixed). The general picture is quite similar to that of the previous section where 
preferential concentration was studied and we see that memory can be important in 
determining the collision rate of inertial particles. 

An interesting difference to the case of preferential concentration is that for large g 
the influence of memory on collision rates is stronger than on preferential concentration 
(compare figures l8l and ITOl) . For example at p = 1000 and St = 1 there is a very small 
change in D 2 due to memory but a more noticeable one in F (8% relative change). 
The collision rate in turbulence is generally believed to be influenced by two effects: 
prefer ential concentration and caustics (also called sling effect), see e.g. Vofikuhle et al\ 
(201J). Therefore this larger increase in the collision rate, as compared to preferential 
concentration, might be due to an influence of the history force on caustics; or more 
specifically, on the relative velocity of colliding particles. 

As detailed in the last section, we are not able to discern the roles of r and g when 
St is fixed, as is the case in figure [TO] This can be only done by looking at the change 
of r as a function of r and g, as has been done for the correlation dimension in figure [Q] 
In the analogous figure for rmemory/Fno memory (not sliown here) one observes a similar 
picture: the influence of memory depends on both r and g. The summary of the findings 
for collision rates is similar as well: for a fixed Stokes number the influence of memory 
increases monotonically when the particles density decreases or (equivalently) when the 
particle size increases. 
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Figure 10. The collision kernel F as a function of t he Sto kes number for different densities with 
and without memory. It is normalize d by Fp = ^/87^/15 ( 2r)^/T n the collision kernel based on 
the dynamics of tracers, as derived bv ISaffman fc Turner! lll956fl . 


8. Comparison with the work of lOlivieri et al.\ (j2014h 

Recently lOlivieri et al. (2o3) published a study of the history force in turbulence. 
During the preparation of the pre sent work discrepanc ies with this publication have been 
found. The type of flow used by Olivieri et al\ ( 2014 1 is very similar to ours, therefore 
we present here a detailed quantitative comparison. We note that there is a difference 
in the methods of forcing the turbulence, which might be a cause of the discrepancies. 
Our forci ng (which is d e scribe d in appendix |D|) varies on a large time scale while the 
forcing bv lOlivieri et all (20l3) changes randomly at every time step (L. Brandt, private 
communication, 2014), i.e. the flow and the particles are forced on small time scales. 

Olivieri et al. ( 2014ll analyse the contributions of different forces to the whole acceler¬ 
ation of the particle. In our notation their quantities are 
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ail ax ai/ux 


Figure 11. PDFs of the ratios Ui^xjax, i € {PG, SD, BA} for g = W and St' = 0.1 without (a) 
and with (b) memory. The vertical dashed lines show the positions of the averages of the PDFs. 
The average values predicted by the approximate formulas (I8.10ll - (l8.12ll are shown with vertical 
dotted lines. 
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Tp 



v) 


( 8 . 1 ) 

( 8 . 2 ) 

(8.3) 

(8.4) 


representing the pressure gradient, the added mass, the Stokes drag and the (Basset) 
history force. The factor ^ = Tp/Tp = (2 p+l)/(2p) ac counts for the slightly different 
definition of the particle response time bv lOlivieri et all ( 2014 1: = 2gr'^f{9v). We will 

use the corresponding Stokes number St' = Tp/Tp = St/^ for the particle parameters in 
this section. Note that f and /3 are shorthand notations for certain combinations of g and 
that these three quantities are equivalent to each othen _ 


To quantify the contributions of the different forces lOlivieri et all (|2014l l consider the 
component-wise ratios of the accelerations given above and the full particle acceleration 
a = dv/dt = Mpc -b Oam + 0,50 + Oba: 

(8.5) 




Because of the statistical isotropy of the underlying turbulent flow, it suffices to consider 
only one component (we choose the x-component). We omit here the contribution of the 
added mass effect as it is fully determined by pressure gradient contribution: aAM,a;/«a; = 
2 {OpG^x/Ox T 1/i?)- 

Figure [TT] shows the PDFs of the three ratios (18.51) for g = 10, St' = 0.1 in our 
simulations with and without memory. (Note that the choice of c olours in the figures of 
this section is different from rest of the paper; it is the same as in ( Olivieri et al. (20l3)), 
to facilitate the comparison). Let us consider the case without memory first. It is shown 
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much smaller. Note that the differences appear also for the case without memory and 
thus cannot stem solely from the treatment of the history force. 

Let us get an independent estimate of the ratios (|8.5|l . To this end we use the approx¬ 
imation m and the weakly inertial estimate for the slip velocity 


u — V 




( 8 . 6 ) 


(which can be obtained by expandin 
see e.g. lDruzhinin fc Ostrovskvl ( 199 
acceleration 


mi in powers of Tp and keeping the linear terms, 
) and obtain estimates for (18.31) , (18.41) and the total 



(8.7) 

r Di/, 

1-/3 — 
r] Dt 

(8.8) 

( r \ D?i 


(8.9) 


The case without memory is contained for a = 0. For the ratios (18.5|) we get estimates 
which depend on the particle parameters and a only: 




31 +(1-/3) 

ai-/3) 

l + a^(l-/3) 

?«^(l-/3) 

l + aHl-/3) 


( 8 . 10 ) 


( 8 . 11 ) 


( 8 . 12 ) 


These estimates are shown in figure 1111 as dotted lines. For the case without memory 
they coincide perfectly with the mean values obtained from the simulations (the dashed 
and dotted lines overlap). With memory the agreement between the estimate and the 
measured values is not perfect but still good. The main source of error here is probably 
the value of the coefficient a. We t hus see that t h e abo ve estimates fit quite well to our 
simulations but not to the ones bv lOlivieri et all ( 2014 1. 

Let us come to the case of neutrally buoyant particles. Figure [12^ shows the PDFs 
of the ratios (j8.5p for g = 1 and St' — 0.01 with memory (without memory the picture 
is the same). The PDFs are sharp peaks around the values 0 and 1. We find that the 
contribution by the Stokes drag and the history force are negligible and basically the 
full contribution to the acceleration is from the pressure gradient Opc. This is expected 
for neutrally buoyant particles with very small Stokes numbers (St' = 0.01 here) as they 
should behave like tracers: with v = u we get = ciba = 0, a = Opc = D w/Df and thus 
0'STi,xl<ix = o-BA^x/cix = 0, OpG.x/fla; = 1- Tlus picture is very different in ( Olivieri et all 
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Figure 12. (a) PDFs of the ratios ai^^la^ i € {PG, SD, BA} for g = 1 and St' = 0.01 with 
me mory, (b) The ratio H8.13II as a function of the particle size from our simulations and those 
bv lOlivieri et al\ ll2014ll . 


||20141 }. see their figure 4. The average contribution of the Stokes drag is far from zero 
and even larger than the contribution of the pressure gradient. Another strong difference 
is that the PDFs are rather broad in contrast to the sharp peaks found here. 

To obtain an overview of the deviations (for p ^ 1) let us consider the ratio 


{^BA,x/^x} 
{^SD,X /^x) 


(8.13) 


with the aim to describe the importance of the history force relative to the Stokes drag, 
akin to (an) / (as) studied in section |3l Figure [12b shows the ratio (j8.13ll obtained from 
our simulations (the dependence on r/ p is similar to t h at in figure [2b) and compares these 
with three data points obtained from ( Olivieri et all ( 20141 1'). We see that the difference 
is considerable, the largest deviation is by almost one order of magnitude (for g = 10, 
St' = 0.1 and thus r/rj = 0.21). 

To quantify preferential concentration Olivieri et al. { 2o3) use the radial distribution 
function g{R) and find that memory reduces g{R) and thus preferential concentration. 
Generally this is what we also find, but the amount of reduction is much smaller in our 
case. For example, at St' = 1 we find the relative reduction 


1 - 


g{R = r/p) 


memory 


g{R = r/p) 


(8.14) 


no memory 


to be 25% for p = 10 and 3.5% for p = 1000 while [Olivieri et all (120141 ) find 52% and 
23%, respectively. 

To conclude , we se e that there are strong differences between our results and those by 
Olivieri et all ( 2014l l. A comparison with analy tical estimates of th e forces favours our 
results. From our perspective the simulations bv lOlivieri et al. ( 2014 1 seem to overestim¬ 
ate thejniDortance_oftheJiistor^Torce^^WjioteJioweverJhal^he_conclusiona^^ Olivieri 
et al. (l2014h are similar to ours on a more general, qualitative level. 


9. Summary and discussion 

In the present paper we we have studied the importance of the the history force - a 
memory effect - for the motion of inertial particles in a turbulent flow. We have found 
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that it can be quite important (depending on the particle parameters). Its effect is to 
reduce (i) the slip velocity, (ii) the difference in acceleration between inertial particles and 
tracers, (iii) the preferential concentration and (iv) the collision rate. We can concisely 
summarize these finds as follows: the history force causes inertial particles to stay closer 
to the flow and to behave more like trace rs. This is in accordance with previous findin gs 
in a smooth flow with chaotic advection ( Daitche fc Tel ( 201lh : Daitche fc Tell ( 2014l l'). 


How large the effects of memory are depends on the particle parameters. We have 
found simple approximative relations for the magnitude of the history force relative to 
the Stokes drag 


^history force qJ" 
^Stokes drag ^ 


and for the reduction of the slip velocity by memory 


Ju 

|u- 


- 7i 

I memory 

11 1 

I no memory 


1 

l + a^' 

V 


(9.1) 


(9.2) 


The constant a is expected to be on the order of unity and from numerical simulations 
we have determined a — 0.69. These relations make it possible to quickly obtain an 
estimate of the importance of the history force (when a is not known, setting a = 1 
should still yield a reasonable estimate). The fact that these relations fit well to numerical 
simulations shows that for the magnitude of the history force and for the reduction of 
the slip velocity the determining parameter is the particle size and that the density is of 
minor importance. 

This, however, is not true for all effects of memory. We have found the influence of 
the history force on particle acceleration, preferential concentration and collision rates to 
depend on both the particle size and density. In these cases we were not able to find simple 
estimates as above. A general summary of the findings on preferential concentration and 
collision rates is that for a fixed Stokes number a decreasing density or equivalently an 
increasing particle size lead to stronger effects of the history force. These effects can 
become quite strong for small densities/large sizes. 

An important field of study for inertial particles is the motion of water droplets in air 
(e.g. in a cloud), where the density ratio g r; 1000 is very large. In this case we have 
found the effect of memory to be rather small. This is also true for the relative magnitude 
of the the history force and the reduction of the slip velocity, which do not depend on 
the particle’s density (at least approximately). At first this might seem contradictory 
but becomes clear when we note the following: we considered a fixed range of Stokes 
numbers {St G [0.1,3]), the corresponding particle sizes for g = 1000 are quite small 
{r/f] G [0.02,0.12]) and hence also the effect of memory. It might, however, become more 
significant for larger r/r]. Even within the investigated range, an interesting hnding is 
that the change of the collision rate caused by the history force at p = 1000 and = 1 is 
8 %. Although this is still small, it is significantly larger than the effect on all other studied 
quantities for g = 1000. It remains to be seen whether this has any practical implications. 


Valuable discussions with T. Tel, M. Wilczek, L. Brandt and J. Bee are acknowledged, 
as well as the support by the Studienstiftung des deutschen Volkes. The Eulerian part of 
the simulation code has been written by M. Wilczek. 
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Appendix A. Faxen corrections in comparison to the history force 

The aim of the present section is to compare the magnitudes of the Faxen corrections 
and the history force. We will see that the former is much smaller than the latter for 
most particle parameters considered in this paper. 

With the Faxen corrections the evolution equation becomes 


dt ^Dt 


!1au]-^ v-u-Cau]-M $-] 


10 J Tp 




V — u - Am 

6 


_ (Al) 

Gatignol 1 1983h : Maxev &: RilevI (198^. The full Faxen corrections are actually 


see 

defined in terms of surface and volume integrals of u and Da/Dt. The Am terms in (lA II) . 
which are the commonly used forms of Faxen corrections, are the first non-vanishing terms 
of a Taylor expansion of these integrals. We restric t ourself here to these ap proximations 
(the full Faxen corrections have been considered bv ICalzavarini et al\ (2012 h in the case of 
neutrally buoyant particles through a special numerical scheme). The second limitation of 
our comparison is that we do not solve (lA II) but equation (11.21) without Faxen corrections 
and only afterwards evaluate Am along the particle trajectories. Our aim is to obtain 
a rough comparison of the typical magnitudes of the Faxen corrections and the history 
force; we expect that our approach, in spite of these limitations, is sufficient for this 
purpose. 

In section |3] we have compared the history force to the drag and found a simple approx¬ 
imate rule for their ratio (j3.2l) . We proceed in the same way for the Faxen corrections, 
by considering the ratio to the Stokes drag. The last two terms in yield the ratio 

^ |Am| / \v — m|. By approximating the slip velocity with 


V — u 

which follows from (14.11) with ^ r 

|Am| 

6 |m — m| 


_ 1 Dm 

1 + Dt ’ 

V 


Du 
Dt ’ 


and using (jl.4ll we obtain 


1 13 


2\p-l\ 


1 -I- a- 


r\ |i^Am| 


r]J |DM/Dt| 


(A 2) 


(A3) 


Approximation (lA 21) is similar to ()8.6I) . but is of the slightly higher order 3/2 in Tp (as 
1/ (l-|-ar/ry) r; I — arjrj and r oc y/ff). By averaging over tracers in our simulation we 


find 


(I^Am|) 

(|DM/Dt|) 


0.14. 


(A 4) 


Because an average over tracers is equivalent to an Eulerian average, this ratio can be 
considered a property of the underlying flow, independent of the particle picture. Also, 
this results implies that the main contribution to Tdu/Ot = — Vp -I- ^Am is from Vp and 
thus justifies the commonly used name “pressure gradient” for Ou/Ot. Figure [13] shows 
the ratio ^ (| Am|) / (|m — m|) obtained from simulations together with the approximation 
(|A 3|) for the numerical value (Eli). It fits quite well for g < 10. For p > 10 there is a 
notable mismatch, nevertheless (IA3I) provides a reasonable order of magnitude estimate. 

Figure US! also shows the ratio of the history force and the Stokes drag {an) / (os), 
which was studied in sectionjS] We see that for all studied parameter values the magnitude 
of the Faxen corrections is smaller than the magnitude of the history force; in most 
cases it is much smaller. From all the considered densities the relative size of the Faxen 
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Figure 13. A comparison of the magnitudes of the Fa:xen corrections and the history force. 
Shown are the ratio ^ (|Au|) / (|v — n|) (filled symbols), its approximation (lA311 (dotted lines) 
and the ratio {an) / {as) (unfilled symbols). 


corrections are largest at p = 0.5. This is expected because this value is closest to g = 1, 
in which case v k, u and thus the ratio ^ |Am| / \v — u\ becomes very large. 

For for p = 1 we can expect an ~ 0 (see also appendix |B|) and because the Faxen 
corrections remain finite, they are expected to be much more important than the history 
force for neutrally buoyant particles (indeed even more important than the Stokes drag). 

If we use the approximation (lA 3|) and extrapolate to smaller particles sizes, we see 
that the Faxen corrections become larger than the history force when r is small enough. 
Indeed the ratio ^ |Ait| / \v — u\ tends to a constant for r —>■ 0 while an/as goes to 
zero, see figure 1131 Thus for small r (how small, depends on g) the Faxen corrections 
become more important than the history force. However this is not the case for all of our 
studied parameter values. Also, when g is not close to 1 both terms are quite small when 
the Faxen terms become more important; then both might be neglected. 

Let us now come to the Faxen correction from the first term in (lA II) . Using (lA 21) we 
find the approximation 


^\v-u\^ lO\P-l\\r]J ^ r]J |DM/Dt| 

' P 


(A 5) 


and for the unknown ratio 


= 0.11 

(|Dti/Dt|) 


(A 6) 


from simulations. Using this value and a = 0.69 (obtained in section |3]) we can estimate 
that the Faxen correction from the first term in (lA II) is smaller than that form the second 
or third term when r/rj < 2.5. Furthermore we find that the ratio (lA 51) is smaller than 
{an) / (as) for all the considered parameter values (i.e. those shown in figure [T5]) . 

We thus conclude that the magnitude of the Faxen corrections is smaller than the 
magnitude of the history force for all the investigated particle parameters; in most cases 
it is much smaller. 
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Appendix B. Note on neutrally buoyant particles 

Up to now we did not consider the case of neutrally buoyant particles. This case is 
special because trajectories of tracers, i.e. v — u, are valid solutions of (11.21) for g = 1. 
Indeed we have found in our simulations (for Stokes numbers 0.01, 0.1 and 1) that neut¬ 
rally buoyant particles behave very similar to tracers: the mean slip velocity is practically 
zero, the acceleration PDF is the same as for tracers, there is no preferential concentra¬ 
tion and collision rates equal that of tracers (with an artificially assigned, corresponding 
radius). The inclusion of the history force does not change any of these findings and thus 
seems not to be important for neutrally buoyant particles. 

We note that for larger particle si zes (not considered here) neutrally buoyant particles 
can behave differently from tracers; Calzavarini et all ( 2012l l have studied the effect of 
memory in this case (along with the influence of Faxen corrections and non-linear drag). 

Concerning the forces acting on a particle, we find fag) ~ 0 and («//) ~ 0 for all 
studied Stokes numbers (see also figure [T^K). Thus, the overwhelming contribution to 
the particle acceleration comes from the pressure gradient ap. This does not mean that 
the Stokes drag is unimportant, because neglecting it (and the history force) in (11.21) 
would allow for solutions of the type v = u + c with an arbitrary constant c, which 
is undesirable. One needs the Stokes drag as a restoring effect which keeps the particle 
close to a tracer trajectory (even though it has to act rarely). The history force can be 
also considered as a restoring effect (indeed it is sometimes called the unsteady drag). 
The fact that we do not find any influence of memory suggests that the Stokes drag is 
already sufficient to keep the particles very close to the tracer trajectories. 

Note that we neglected the Faxen corrections in our simulations. This might be a 
serious limitation for neutrally buoyant particles as the findings of appendix El suggest. 
Thus the results of this section should to be treated with some caution. 

A fundamental question c oncerning neut r ally b u oyant particles is whet her the tracer 
solutions are stable, see e.g. Babiano et al. ( 200(11 1: Sapsis & Haller ( 2008ll. The history 
force m ight play a very non-trivial role here. For example it was shown bv IPaitche fc Tel 
(201 ill that the history force can destroy attractors and that the convergence towards 
attractors becomes very slow (^ with memory. Thus a thorough study of the 

influence of memory on neutrally buoyant particles would have to include a stability 
analysis of tracer trajectories; this, however, is beyond the scope of this paper. 


Appendix C. A time-stepping scheme for particles with memory 

Let us start by rewriting the evolution equation (11.21) for the slip velocity w = v — u: 


dir , ^ , da ^ „ 1 

= Pw-Vu - w - 


dt 


dt 



dx 

dt 


— w + u, 


(Cl) 

(C2) 


where du/dt = dtu + v ■ Vu. We wish to obtain a numerical solution at the time points 
tn = nh, where h is the time step. Integrating (1C II) from to tn+i we get 


Wn+l = Wn + - 1) {Un+l - Un 



dr 


(C3) 
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where ^ = —y^3/3 /(ttTp). To obtain a numerical stepping scheme we have to approximate 
the integrals. We begin with the first two: 


rtn + 1 ^ 

/ /(r)dT Ri h'^xf^''f{Tn+i-i) 

i =0 


(C4) 


The coefficients are well known. For an explicit scheme, i.e. Aq = 0, they are the 
coefficients of the Adams-Bashforth schemes: 

A ^’2 = ( 0 , 1 ) 


( 0 ,^,- 


3 1 


= ( 0 , 


2 ’ 2 , 

23 4 5 


(C5) 

(C6) 

(C7) 


For an implicit scheme they are the coefficients A^^’™ of the Adams-Moulton schemes: 


xM,2 _ / 1 1 

“ 2 ’ 2 


5 8 1 


A"”" = —,- 

' 12 ’ 12 ’ 12 


(C8) 

(C9) 


The errors of these approximations are for A'®’"* and for A'^’™. The 

history integrals in l|C 3|) are approximated with 

rtn " 

zdr Ki y/h'^ (CIO) 




i=0 


where = w{tn) and the coefficients are given by lDaitch3 ( 2013h . We thus obtain 
the following scheme for the slip velocity: 

m 7 m 

Wn-\-l = Wn-\- i/3 - 1) (Wn +1 “ Af [w ■ ■-^ 


2=0 


2=0 


/ 22+1 


E ^fiiWn-i\ (Cll) 


\2 = 0 


2=0 


Note that the right-hand side contains Wn+i, but only linearly. Bringing it to the left- 
hand side yields an explicit scheme: 


^1 + Wn+l = Wn + - 1) {Un+l - U„) 


-Ph - — E 


M.m 


72 + 1—2 


2 = 0 


2=1 


( 72 + 1 72 A 

E /^2'^72+l—2 ^ ^ 


2=0 


^72+1 - ^72 h ^ ^ A.^ iWn-\-l — i “h 


(C12) 


(C13) 


2=0 


Using Af’^, Af^’^ and the third-order coefficients pi from ( Daitchel ( 2013h l the above 
scheme has a one-step error of O(h^), i.e. it is a third-order scheme. 
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As this scheme relies on previous values of x and w, one needs to start the integration 
with lower order schemes: the first step with Afand the first-order ptf, the second 
step with A^’^ and the second-order 

To further improve the accuracy of the scheme, corrector steps have been applied. A 
corrector step is performed by predicting a;„+i, Wn+i with the explicit scheme (1C 121) - 
(|C 13|) . replacing the Xf in this scheme with Af^ (thus making it implicit), inserting the 
predicted Wn+i in the right-hand side and calculating a corrected estimate for 

Xn+i, lUri+i- (The corrector steps are not essential for the numerical scheme.) 

That the one-step error of this scheme scales indeed with h^, has been verified in a case 
where an analytical solution o f (jl.2l) is av ailable: a partic le moving in a two dimensional 
vortex ( Candelier et al. (2004)), see also ( Daitche (2013)). 


Appendix D. Forcing the turbulence 

We solve the vorticity equation in a triply periodic box of length Lbox = 27r. In Fourier 
space the equation becomes 

dtiOk = ik X iF[u X tjj]k — vk^uik + Fk- (D 1) 

Here a; = V x m is the vorticity, fc G the wave vector, F the Fourier transform and 
F the forcing. The latter is implicitly defined as described in the following. We force by 
inserting the energy lost during one time step back into large scale modes ijJk with k in 
the forcing band 

H = {fc |3 < < 9} . (D2) 

Let be the current vorticity and the vorticity after one time step without forcing 
(more precisely, when using a scheme consisting of multiple stages, like a Runge-Kutta 
method, is the vorticity after one stage). The forcing is applied by rescaling the 

modes in the forcing band: 


+1 kGB 

" k^B 

(D3) 

The factor / is chosen such that the energy of the flow 



(D4) 

remains constant, i.e. 


E[u:l+^]=E[u;l]. 

(D5) 


The modes{a;fc| k ^ B} evolve freely and the energy is kept constant by inserting energy 
into the modes {wfej k G B}. Note also that only the modulus of the modes is modified, 
their direction can evolve freely, which improves the statistical isotropy of the flow. 
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